function [loglike,TCPmodel,EQD2,Effect,ps] = loglikelihood_poisson(par_long,data_TCPexp,inds,N,dep_a,dep_b,treatment_type)
%This function calculates the fitting loglikelihood of a Poisson TCP model using some input parameters (par_long and N) and
%a SF expresssion with functional dependencies dep_a and dep_b for every cohort of a clinical treatment database (data_TCPexp).

%Input:

%        1. par_long: TCP model initial parameter values par_long = [alpha alphabeta a b g lambda_prolif t_kickoff TCPmax]
%        2. data_TCPexp: database array with information of clinical RT treatment for several patient cohorts 
%        3. inds: Cohorts in the database with more than 170 pts.
%        4. N: Number of initial tumor cells
%        5. dep_a: functional dependence with dose of the alpha term in the expression of the surviving fraction
%        6. dep_b: functional dependence with dose of the beta term in the expression of the surviving fraction
%        7. treatment_type: String with treatment site ('brain' for brain metastases and 'lung_' for NSCLC.



%Output: 
%        1. loglike: loglikelihood of the best fit
%        2. TCPmodel: TCP model values associated to each cohort of the database calculated using dep_a and dep_b model functional dependencies and the input parameter values         
%        3. EQD2: EQD associated to each cohort of the database calculated using dep_a and dep_b model functional dependencies and the input parameter values
%        7. Effect_opt: Effect,log(SF), associated to each cohort of the database, calculated using dep_a and dep_b model functional dependencies and the input parameter values         
%        4. ps: fitting probability (use of binomial distribution) 


%Model parameters:
alpha=par_long(1);
alphabeta=par_long(2);
a=par_long(3);
b=par_long(4);
g=par_long(5);
lambda=par_long(6);
tk=par_long(7);
TCPmax=par_long(8);
%Treatment information from every cohort in the clinical database: number of treatment schedules, number of fractions, dose per fraction, number of patients/metastases involved, treatment time and tumor control probability. 
schedules=data_TCPexp(1,7); %number of treatment schedules
n=data_TCPexp(:,7+1:2:7+(2*schedules-1)); %number of fractions
d=data_TCPexp(:,7+2:2:7+2*schedules); %average dose per fratcion
pacs=data_TCPexp(:,1); %number of treated patients for NSCLC/ treated metastates for brain
pTCPexp=data_TCPexp(:,3); %Agresti-Coul probability associated to experimental Kaplan TCP value
t_trat=data_TCPexp(:,6); %Treatment time for calculation of accelerated repopulation effect
num_datapoints=length(data_TCPexp(:,7));
vector_ones=ones(max(data_TCPexp(:,1)),1);

%Calculation of treatment effects and EQD2 values.
Effect_central=sum(n.*d.*alpha.*(1+a.*d.^dep_a+d.*(1+b.*d.^dep_b)/alphabeta),2)-lambda*max(0,t_trat-tk);
    if dep_a==-1 %Treatment effect under LQL model
        Effect_central=sum(n.*alpha.*(d+2.*(a.*d+exp(-a.*d)-1)./(alphabeta.*a^2)),2)-lambda*max(0,t_trat-tk);  
    end  
[EQD2]=EQD2_calculator(Effect_central,[alphabeta a b 0 0 lambda tk alpha],dep_a,dep_b);  

%if EQD2 have no NaN values, calculate: TCP, fit probability and  fit loglike (lines 189 to 257) 
    if any(isnan(EQD2))==false 
         muestra=1000;
         %Definition of alpha radiosensitivity values array (alpha_vect) to account for interpatient hetereogeneity
         alphamin_i=max(0,alpha*(1-5/sqrt(g+1)));
         alphamax_i=alpha*(1+5/sqrt(g+1));      
         alpha_vect=[alphamin_i:(alphamax_i-alphamin_i)/muestra:alphamax_i];
         %Definition of array with gamma pdf distribution values associated to alpha_vect values
         pdf_vect=gampdf(alpha_vect,g+1,alpha/(g+1));
         %Calculate treatment effect for alpha_vect sensitivity values for every cohort
               for i=1:num_datapoints
                   Effect1=max(n(i,1).*d(i,1).*alpha_vect.*(1+a.*d(i,1).^dep_a+d(i,1).*(1+b.*d(i,1).^dep_b)/alphabeta)-lambda*max(0,t_trat(i)-tk),0);
                                if treatment_type=='brain'
                                    Effect2=max(n(i,2).*d(i,2).*alpha_vect.*(1+a.*d(i,2).^dep_a+d(i,2).*(1+b.*d(i,2).^dep_b)/alphabeta)-lambda*max(0,t_trat(i)-tk),0);                   
                                    Effect_t(:,i)=Effect1+Effect2; 
                                else
                                    Effect_t(:,i)=Effect1; 
                                end
                      if dep_a==-1 %LQL model
                        Effect1=max(n(i,1).*alpha_vect.*(d(i,1)+2.*(a.*d(i,1)+exp(-a.*d(i,1))-1)./(alphabeta.*a^2))-lambda*max(0,t_trat(i)-tk),0);  
                               if treatment_type=='brain'
                                   Effect2=max(n(i,2).*alpha_vect.*(d(i,2)+2.*(a.*d(i,2)+exp(-a.*d(i,2))-1)./(alphabeta.*a^2))-lambda*max(0,t_trat(i)-tk),0);  
                                   Effect_t(:,i)=Effect1+Effect2; 
                               else
                                    Effect_t(:,i)=Effect1; 
                               end
                     end
                  SF(:,i)=exp(-Effect_t(:,i)); 
               end
         anchura=alpha_vect(2)-alpha_vect(1); 
         %Numerical integration of TCP and Effect considering interpatient heterogeneity
         TCPmodel=TCPmax.*sum(pdf_vect'.*exp(-N.*SF))'.*anchura; 
         Effect=sum(pdf_vect.*Effect_t',2).*anchura;
         %Calculation of the probability that the experimental TCP value yielded the number of controls estimated from the radiobiological
         %model TCPmodel. Use of binomial probability, where factor is all x (tumor control cases)-combinations of a set Ns patients or metastases.
         factor=gamma(pacs+1)./(gamma(pacs.*TCPmodel+1).*(gamma(pacs-pacs.*TCPmodel+1)));
         %if controls=0, factor=1
         inds_with_0controls=inds(round(TCPmodel(inds).*pacs(inds))==0);
         factor(inds_with_0controls)=1;
         %delete inds for those patients
         idx=ismember(inds,inds_with_0controls);
         inds(idx)=[];
         %if control=1, factor=patients
         inds_with_1control=inds(round(TCPmodel(inds).*pacs(inds))==1);
         factor(inds_with_1control)=pacs(inds_with_1control);
         idx=ismember(inds,inds_with_1control);
         inds(idx)=[];
         %binomial fit probability calculation  
         ps=factor.*pTCPexp.^(pacs.*TCPmodel).*(1-pTCPexp).^(pacs.*(1-TCPmodel));
         
         %To avoid divergence (Inf) when patient/metastases number (Ns) is very large, round x (with small effect) and calculate 
         %factor=n!/(x!(n-x)!)=prod([max(x,n-x)+1:1:n))/factorial(min(x,n-x)) with a=max(x,n-x)and b=round(min(x,n-x)); 
             for ij=1:length(inds) 
                   xr=pacs(inds(ij)).*TCPmodel(inds(ij));
                   x=fix(xr);  
                   %If any of the values goes to infinite, use binopdf
                   % first interpolation point (x-1), avoiding negative values
                   x0=max(0,x-1);
                   pacs_i=pacs(inds(ij));
                   pTCPexp_i=pTCPexp(inds(ij));
                   n_mx=pacs_i-x0;
                   a_f=max(x0,n_mx)+1;
                   b_f=min(x0,n_mx);
                   vector_j=a_f:1:pacs_i;
                   vector_num=vector_ones;
                   vector_num(1:length(vector_j))=vector_j;
                   vector_den=vector_ones;
                   vector_den(1:b_f)=1:1:b_f;
                   factor_1=prod(vector_num./vector_den);
                   p1=factor_1.*pTCPexp_i.^(x0).*(1-pTCPexp_i).^(n_mx);

                   if p1<Inf
                   %second interpolation point (x), central point (c)
                      n_mx=pacs_i-x;
                      a_f=max(x,n_mx)+1;
                      b_f=min(x,n_mx);
                      vector_j=a_f:1:pacs_i;
                      vector_num=vector_ones;
                      vector_num(1:length(vector_j))=vector_j;
                      vector_den=vector_ones;
                      vector_den(1:b_f)=1:1:b_f;
                      factor_c=prod(vector_num./vector_den); 
                      pc=factor_c.*pTCPexp_i.^(x).*(1-pTCPexp_i).^(n_mx);

                      if pc<Inf
                         % third interpolation point (x+1), avoiding x+1>pacs 
                         x1=min(pacs_i,x+1);
                         n_mx=pacs_i-x1;
                         a_f=max(x1,n_mx)+1;
                         b_f=min(x1,n_mx);
                         vector_j=a_f:1:pacs_i;
                         vector_num=vector_ones;
                         vector_num(1:length(vector_j))=vector_j;
                         vector_den=vector_ones;
                         vector_den(1:b_f)=1:1:b_f;
                         factor_2=prod(vector_num./vector_den);
                         p2=factor_2.*pTCPexp_i.^(x1).*(1-pTCPexp_i).^(n_mx);

                         if p2< Inf
                             p_v=[p1 pc p2];
                             %Avoid repeated values
                             [x_unique,i_unique]=unique([x0 x x1]);
                             p_unique=p_v(i_unique);
                             ps(inds(ij))=spline(x_unique,p_unique,xr);
                              if ps(inds(ij))< 0
                                 ps(inds(ij))=spline([x x1],[pc p2],xr);
                              end
                         else
                          ps(inds(ij))=calc_binopdf(xr,pacs_i,pTCPexp_i);
                         end
                      else 
                       ps(inds(ij))=calc_binopdf(xr,pacs_i,pTCPexp_i);
                      end
                   else 
                     ps(inds(ij))=calc_binopdf(xr,pacs_i,pTCPexp_i);
                   end  
             end
         loglike=-log(prod(ps));
    else %If there is no EQD2, discard parameters (loglike set to large value, and other variables = NaN or 0)
         TCPmodel=NaN(num_datapoints,1);
         EQD2=NaN(num_datapoints,1);
         loglike=1e100;
         ps=zeros(num_datapoints,1);
         Effect=NaN(num_datapoints,1);
    end
end

function [p]=calc_binopdf(xr,n,p)
x=fix(xr);
x0=max(0,x-1);
x1=min(n,x+1);
p1=binopdf(x0,n,p);
pc=binopdf(x,n,p);
p2=binopdf(x1,n,p);
p_v=[p1 pc p2];
%Avoid repeated values
[x_unique,i_unique]=unique([x0 x x1]);
p_unique=p_v(i_unique);
p=spline(x_unique,p_unique,xr);
%Avoid interpolation of a negative value
if p< 0
    p=spline([x x1],[pc p2],xr);
end
end